Large-scale mRNA transfer between Haloxylon ammodendron (Chenopodiaceae) and herbaceous root holoparasite Cistanche deserticola (Orobanchaceae)

Summary Exchanges of mRNA were shown between host and stem parasites but not root parasites. Cistanche deserticola (Orobanchaceae) is a holoparasitic herb which parasitizes on the roots of woody plant Haloxylon ammodendron (Chenopodiaceae). We used transcriptome sequencing and bioinformatic analyses to identify nearly ten thousand mobile mRNAs. Transcript abundance appears to be a driving force for transfer event and mRNA exchanges occur through haustorial junction. Mobility of selected mRNAs was confirmed in situ and in sunflower-Orobanche cumana heterologous parasitic system. Four C. deserticola →H. ammodendron mobile mRNAs appear to facilitate haustorium development. Of interest, two mobile mRNAs of putative resistance genes CdNLR1 and CdNLR2 cause root-specific hypersensitive response and retard parasite development, which might contribute to parasitic equilibrium. The present study provides evidence for the large-scale mRNA transfer event between a woody host and a root parasite, and demonstrates the functional relevance of six C. deserticola genes in host-parasite interactions.


INTRODUCTION
Increasing evidence suggests the importance of macromolecules such as proteins and mRNAs in both short-and long-distance communication in plants. 1 Plasmodesmata are considered as short-distance transport channels, whereas vascular system carries out long-distance molecular transport between distal tissues. 2,3 mRNA was reported to transfer in phloem. 4 Recent studies have shown that the stem-loop structure of tRNA origin endowed mRNA with long-distance transport capacity. 5 AtRRP44A, a subunit of the RNA exosome, was reported to interact with plasmodesmata and to mediate the cell-to-cell trafficking of KNOTTED1 (KN1) mRNA. 6 Parasitic plants account for about 1% of flowering plants. The most important feature of parasitic plants is a specialized structure called haustorium, which establishes physical and physiological connection channels between the host and parasitic plants and so dominates most of their interactions. 7 Parasitic plants rely on the host to sustain their growth, absorbing water and the nutrients such as photosynthates, amino acids and other intermediate metabolites from the host through the haustorium. 8 Haustoria also function as the connection bridge between the hosts to transmit signaling molecule such as herbivore signal. 9 It was also shown that the biomolecules including proteins and RNAs exchange bidirectionally through host-parasite connections. The earliest example of RNA transfer between the host and parasitic plants is the transmission of RNA virus from the infected host to the uninfected plant through dodder. 10 Small RNAs (sRNAs) such as small interference RNA (siRNA) and microRNA (miRNA) were shown to move between host-parasite to regulate the target gene expression in recipient organisms. [11][12][13] As for mRNA transfer, mobile mRNAs from the hosts tomato, pumpkin and alfalfa to the parasite Cuscuta chinensis have been demonstrated for several genes in the initial test. 14 Through microarray analysis, Roney et al. found that 474 mRNAs were transferred from tomato to dodder. 14 However, large-scale identification of transfer mRNAs was achieved only on the uses of next-generation sequencing technology. Kim et al. used transcriptome sequencing technology to identify large-scale and bidirectional mRNA transfer between the parasite Cuscuta pentagona and the hosts such as Arabidopsis thaliana and tomato. 15 But the functional relevance of transfer mRNAs in host-parasite interactions was not clear.
Orobanchaceae is the largest parasitic angiosperm family in which many are facultative or obligate root parasites. 16 Cistanche is a worldwide genus of holoparasitic desert plants in Orobanchaceae. C. deserticola is a holoparasite which parasitizes on the roots of a woody psammophyte, Haloxylon ammodendron (Chenopodiaceae). 17 It was shown that chloroplast rpoC2 gene was transferred from H. ammodendron to C. deserticola by horizontal gene transfer. 18 Here, we identified nearly ten thousand mobile mRNAs between C. deserticola and H. ammodendron through next-generation transcriptome sequencing and bioinformatic analysis. The mobile mRNAs were ascertained by multiple sequence alignment, PCR validation, and especially utilization of a sunflower-Orobanche cumana parasitic system. The function of mobile mRNAs was also demonstrated for several genes by this heterologous parasitic system. To date, the reports on functional mRNAs exchange between host and parasite were very limited. Therefore, our study provides new insights into the parasitic mechanism from the view point of mobile mRNAs.

Identification of mobile mRNAs between C. deserticola and H. ammodendron
Tight physical connection of host H. ammodendron root and parasite C. deserticola makes it impossible to separate them at the haustorium. Therefore, we combined high-throughput RNA-sequencing (RNA-seq) and stepwise bioinformatic classification to identify the mobile transcripts between host and parasite. Four different set of samples were collected for RNA-seq analyses on Illumina HiSeq2500 platform. These include the root tissue of parasite-free host H. ammodendron (HA), succulent stem of C. deserticola (CD), the root tissue of host H. ammodendron parasitized with C. deserticola (HC), and the haustorial interface (HI). For HC and CD samples, the respective host root and parasite stem were collected 1 cm away from the parasitic junction ( Figure 1A).
Reliable identification of source tissue for a given transcript from interacting organisms is a challenging issue. Therefore Ikeue et al. have developed a useful bioinformatic method to distinguish host and parasite transcripts. 19 Here we used a modified bioinformatic approach which is based on their method to identify the mobile mRNAs between C. deserticola and H. ammodendron ( Figure 1B). A grand total of 723,382,940 clean reads were generated from above libraries (Table S2), and these reads were assembled via different strategies ( Figure 1B). As a result, all ten samples were hybrid assembled into 222,899 unigenes (termed ''Combined'' hereafter, Table S3). C. deserticola and H. ammodendron samples were also assembled into 107,752 and 194,720 unigenes respectively (termed ''Cis'' and ''HAC'' hereafter, Table S3). To visualize the variation as well as the similarity for all samples, we performed a principal component analysis (PCA) on the normalized FPKM values of all the detected genes in ''Combined unigenes''. The PCA plot showed that the data for three biological replicates were clustered closely and were separated in different samples, especially between H. ammodendron and C. deserticola ( Figure S1). Meanwhile, ORF prediction showed that 149,825 (67.23%) of the ''Combined'' unigenes were assumed to encode putative proteins ( Figure 1C). The individually assembled unigenes with the expression threshold value of FPKMR0.3 from four different samples were cross-checked with the ORF prediction results (Tables S4-S6). This analysis identified that 69.01-75.90% of the unigenes from the four samples CD, HC, HA and HI had putative ORF ( Figure 1C). Full-length transcriptome sequencing was carried out and the data, CD_FL for C. deserticola and HA_FL for H. ammodendron, were used for assembly error correction for both host H. ammodendron and parasite C. deserticola as both species lack genomic data (Tables S8-S11). When the assigned unigenes from above four samples were filtered with dual BLASTN (E value = 1e À10 ) against Cis, HAC, CD_FL and HA_FL, 94.78-99.00% of them had reliable consensus sequences ( Figure 1C and Table S7). Based on the above analysis, 17,379 (HA-exclusive 14,810 and HA-containing 2,569) common unigenes were finally retrieved from HC and CD on Venn diagram analysis of the unigenes from CD, HC and HA ( Figure 1D). Among them, the former 14,810 unigenes are most probably of parasite CD origin because of their absence in HA, but a proportion of them might also represent HC/CD mobile mRNAs which are present only in HC but not in HA owing to their up-regulation on parasite attachment. The latter 2,569 unigenes most probably represent unidirectional HC/CD transfer mRNAs as they are also present in HA. The other unigenes were excluded from candidate mobile mRNAs as they were not shared between HC and CD. One should also note that 14,810 common unigenes in [HC, CD] outnumbered 166 common unigenes in [HA, CD] by 89 times, highly suggesting the reliability of our analyses as the latter was logically impossible and was caused by the inevitable errors from unigene assembly or bioinformatic analyses ( Figure 1D).
The origin of putative mobile mRNAs was further confirmed via dual BLAST analyses to assign their sequence origin ( Figure 1B and Table S12). Because of the lack of genomic data for both host and parasite, iScience Article we utilized available sequence data for Chenopodiaceae and Orobanchaceae family organisms to which H. ammodendron and C. deserticola belong respectively. Our hypothesis is that mobile mRNAs of H. ammodendron origin most probably have higher homologies with their orthologs from other Chenopodiaceae species than those from Orobanchaceae species, while it is vice versa for mobile mRNA of C. deserticola origin, and the conserved orthologs probably have high similarity with both Orobanchaceae and Chenopodiaceae. This analysis could help us to ascertain the origins of the mobile mRNAs as they could be assigned to family level via homology searching. To prove our hypothesis, phylogenetic tree of thirty-seven species, including six Orobanchaceae and four Chenopodiaceae species, was constructed using OrthoFinder. The far evolutionary distance between Orobanchaceae and Chenopodiaceae species ascertained the reliability of our hypothesis ( Figure 2A and Table S13). Gene loss analysis showed that the transcripts for many orthogroups and photosynthesis related genes were absent in C. deserticola (Figures 2B   Figure S2). These results not only indicated the parasitic property of this species but also the reliability of unigene assembly from our analyses, and also supported our hypothesis to use relative plant species for confirmation of sequence origin. Therefore, the candidate transfer unigenes from above analyses were separated with dual BLASTN (E value = 1e À10 ) against a collection of public available sequence datasets from Orobanchaceae (transcriptome of Triphysaria versicolor, Striga hermonthica and Phelipanche aegyptiaca; EST and mRNA sequences from NCBI) and Chenopodiaceae (next-generation and full-length transcriptome of H. ammodendron; EST and mRNA sequences from NCBI) ( Figure 1B and Table S12). As a result, 7,496 unigenes were confirmed to originate from the parasite C. deserticola, accounting for 9.66% (7,496/77,615) and 13.70% (7,496/54,721) of the unigenes in destination HC samples and source CD, respectively; and 2,370 unigenes were assigned to host H. ammodendron, accounting for 3.38% (2,370/70,119) and 4.15% (2,370/57,091) of the unigenes in source HC and destination CD samples, respectively ( Figures 2D and 2E, Tables S7, S12, and S15). The common 2,931 unigenes were too homologous to be assigned to the exact source organism. These results indicated that nearly ten thousand (7,496 + 2,370 = 9,866) unigenes could transfer between root parasitic plant C. deserticola and woody plant host H. ammodendron. Much more (7,496/2,370 = 3.16-folds) mobile RNAs were of parasite C. deserticola origin than of host H. ammodendron origin, showing the mRNA transfer bias in a parasite/host direction.

General properties of mobile mRNAs
Although mRNA transfer has been documented between stem parasitic plant in Cuscuta species and their hosts, 14,15,20 the transfer mechanism remains poorly understood. It was hypothesized that mRNA transfer might be driven by specific selective determinants, such as abundance, size, featured motif, or specific modification of the mobile mRNAs. 5,15,21,22 Therefore, we asked whether these characteristics of mobile mRNAs could influence the transfer events between root parasitic plant C. deserticola and host H. ammodendron or not. FPKM analyses showed that transcript abundance of the majority of mobile unigenes of C. deserticola origin was higher than that of non-mobile unigenes from same species. However, the transcript abundance of mobile mRNAs of H. ammodendron origin showed divergent patterns of comparable to higher FPKM values than those of non-mobile mRNAs from same species ( Figure 3A). Moreover, mobile mRNAs of CD origin showed higher abundance in source than in destination tissue, where we speculate that the mobile mRNAs might be diluted or degraded after transfer event ( Figure 3B). Compared with this, mobile mRNAs of HC origin showed no uniform correlation of transcript abundances between source and destination tissues ( Figure 3B). These results demonstrated that mRNA transfer could be correlated with unigene abundance, but the mechanistic basis for such correlations remains obscure.
We assume that mobile mRNAs might be present in haustorium as they should translocate through parasitic interface. Most of the mobile unigenes could be found in haustoria unigene dataset, accounting for 99.21% (7,437/7,496) and 80.89% (1,917/2,370) of the transfer unigenes originated from C. deserticola and H. ammodendron, respectively ( Figures 3C and 3D). This result indicated that root parasitic plant and its host exchanged their transcripts through haustorial junction, which serves as mechanical bridge between host and parasite and mediates the transfer event. Further evidence for haustorium-mediated selective mobility of unigenes came from plots of unigene abundance in haustorium versus abundance of the same unigenes in source C. deserticola and H. ammodendron, respectively. The plot of CD/HC mobile mRNAs showed that expression levels of most unigenes in the source CD sample were positively correlated with those in haustorium ( Figure 3E). This result was similar to previous report about mobile mRNA between Arabidopsis and Cuscuta (Kim et al., 2014), indicating that most unigenes followed the same dynamics of movement. The plot of HC/CD mobile mRNAs showed a more-dispersed pattern of transcript iScience Article abundances without obvious correlation between source HC and haustoria ( Figure 3F). These results demonstrated that root parasitic plant C. deserticola and its host H. ammodendron could exchange their transcripts in a haustorium-mediated selective manner, but the dynamics of movement differed between two directions.
KEGG analyses showed the overall functions of the mobile mRNAs ( Figures 4A and 4B). In the KEGG analysis, the terms related to RNA transport and metabolism, such as ''RNA transport'', ''spliceosome'' and ''mRNA surveillance pathway'' were enriched in CD/HC mobile RNAs, whereas none of them were detected in HC/CD mobile RNAs. The enrichment of above pathways is consistent with the more frequent transfer of mRNA in CD/HC direction than that in the opposite direction ( Figure 4A). ''RNA transport'' related mobile mRNAs in parasite C. deserticola might facilitate the large-scale transfer of mRNAs in the CD/HC direction. Furthermore, the presence of ''spliceosome'' and ''mRNA surveillance pathway'' related mobile RNAs from the parasite into host indicates that these mobile mRNAs might be involved in the maturation of parasite-derived mRNAs, or they may modify the mRNAs of host origin to influence iScience Article the recipient plant's physiology. 23 In addition to ''RNA transport'', ''protein export'' pathway was also enriched in CD/HC mobile RNAs ( Figure 4A). More frequent transfer of proteins than RNAs was observed in the association of stem parasite dodder and its host Arabidopsis and soybean. 24 The term ''autophagy'', a process involved in programmed cell death (PCD) and immune response, was enriched in CD/HC mobile RNAs ( Figure 4A). In HC/CD direction, metabolic pathway genes for protein biosynthesis are significantly enriched compared with CD/HC direction. The term ''ribosome'', ''protein processing in endoplasmic reticulum'', ''ribosome biogenesis in eukaryotes'' and ''aminoacyl-tRNA biosynthesis'' are all related to protein biosynthesis and related genes were highly represented in HC/CD mobile RNAs ( Figure 4B).

Experimental confirmation of selected mobile mRNAs
We tried to test the independent transfer event but it was impossible to experimentally validate all of these mobile mRNAs. Therefore, we randomly selected some of the unigenes to verify the transfer event (Tables 1  and S3-S14). Polymerase chain reaction (PCR) assay was carried out to validate the mobile mRNAs and genomic DNAs (gDNAs) of both host and parasite were used as control templates. We reasoned that the PCR product could only be amplified from gDNA of source organism other than gDNA of destination organism for a specific mobile mRNA, as the mobile mRNA can only be transcribed from source gDNA but not from destination gDNA. To avoid the amplicon size differences between gDNA and cDNA, we used the gene fragments with no putative intron for PCR assay. We confirmed ten CD/HC and three HC/CD mobile mRNAs using this assay (Table 1 and Figures S3-S14). We could amplify HC/CD mobile sequences from HA gDNA, HA cDNA, HC cDNA and less abundantly from HI cDNA and CD cDNA while not from CD gDNA. Similarly, we could amplify CD/HC transfer sequence from CD gDNA, CD cDNA, and less abundantly from HI cDNA and HC cDNA while not from HA cDNA and HA gDNA ( Figure 5). The presence of corresponding nucleic acids for mobile mRNAs in source gDNA but not in destination gDNA confirmed the transfer direction. Less abundant transcript level in destination than in source sample indicates that mobile mRNA is less represented in the recipient sample, which is consistent with the sequencing results (Figure 3B). As a negative control for this PCR experiment, we tested two HA non-mobile mRNAs (HabHLH110, HabZIP60) and two CD non-mobile mRNAs (CdbHLH47, CdMYB3), and we could not amplify them from CD cDNA or HC cDNA respectively ( Figure 5). Above data indicated the reliability of the PCR results.

Further confirmation and functional analysis of CD/HC mobile mRNAs
Based on the above results, we aimed to further validate the transfer event by in planta overexpression and to characterize the function of the mobile mRNAs. We mainly studied CD/HC mobile mRNAs for validation as they were more broadly and confidently identified than HC/CD mobile RNAs. H. ammodendron-C. deserticola parasitization system is not applicable for this purpose as the transgenic approaches are not available for both the host and parasite. Therefore, we utilized sunflower-O. cumana parasitic system in which the latter belongs to Orobanchae species as the same with C. deserticola. First, we evaluated the phylogenetic relationship between root parasitic plants C. deserticola and O. cumana, as well as corresponding hosts H. ammodendron and Helianthus annuus via phylogenetic analyses using rbcL (RuBisCO large subunit) sequence. The evolutionary tree showed that C. deserticola had a close phylogenetic relationship with O. cumana, as was also the same for the two hosts. However, the stem parasitic plant Cuscuta australis had a far phylogenetic relationship with both parasitic plants and their hosts, which indicated the reliability of the heterologous system ( Figure S15). Previously, we established a sunflower transient expression system in our lab by seed-soak agroinfiltration (SSA) method. 25 Here we modified this SSA method by using pBI121 vector for overexpression of the genes for mobile mRNAs. Candidate gene was fused with GFP (green fluorescence protein) sequence and driven under constitutive CaMV 35S promoter for in planta expression in pBI121 vector. The recombinant constructs were transiently transformed into host sunflower via SSA method, and the parasite O. cumana was attached to the root of host seedling. Total RNA was isolated from both host and attached parasite, and the mobile transcript gene and GFP mRNA was detected with RT-PCR approach, in which host signal indicates the expression of candidate genes in the host and parasite signal denotes the successful transfer of mRNAs from host to parasite. An A. thaliana well-characterized mobile transcript gene AtTCTP1 (Translationally Controlled Tumor Protein) was fused with GFP iScience Article sequence to serve as positive control for transfer event and GFP alone served as negative control. 26 Consistent amplification of these genes in transformed but not control host indicates the successful expression of them in the host plant ( Figure 6A). PCR products were detectable in positive control AtTCTP1-GFP fusion group but not in negative control GFP group in the parasite ( Figure 6). It indicates that GFP mRNA alone is not mobile, whereas fusion of mobile mRNA of AtTCTP1 endows it the ability to move from host to parasite. RT-PCR analysis with both candidate gene specific and GFP specific primers confirmed the existence of the mRNAs for the candidate genes, although they were less abundantly accumulated in the parasite than in the host ( Figure 6B). The results indicated the successful transfer of these candidate mobile mRNAs from host to parasite. We should note that the mechanism of CD/HA movement and sunflower/O. cumana movement of CD mRNA could be different, and that directionality is also of concern to experimental validation. However, the above data provide some evidence for the mobility of these mRNAs given the present technical shortage of transgenic approach for the study of host-parasite communication.
We further investigated the role of mobile mRNAs on parasitic event. On the overexpression of above-verified mobile mRNAs, twenty-day-old sunflower seedlings were parasitized with O. cumana, and 45th day after parasitization O. cumana was removed from the soil for further analysis on haustorium development. As shown in Figure 7A, the number of haustorium (tubercle) formed on CdAPM, CdATPA, CdPGI and CdRPL22 groups increased by 78.1-110.5% whereas CdNLR1 and CdNLR2 groups decreased by 70.5% and 28.6% respectively, compared with control GFP group. There was no significant difference between five other treatment groups and the control GFP group. Representative parasites formed on sunflower root were shown in Figure 7B.
CdNLR1 and CdNLR2 negatively regulate parasite development CdNLR1 and CdNLR2 encoded resistance (R) proteins with putative NBS-LRR domain, respectively (Table 1). Given the importance of these two genes, we further confirmed the mobility of their mRNAs by detecting GFP signal in earlier haustorial development stage. On overexpression in the host sunflower, fusion of these two   Figure S16). In later growth stage, overexpression of CdNLR1 and CdNLR2 caused visible tubercle necrosis which was not seen in control GFP group ( Figure 8A). When the tubercle was observed at pre-elongation stage, compatible attachment was observed in the control GFP group, in which successful vasculature connection formed between host and parasite. In contrast, incompatible attachment was observed in CdNLR1 and CdNLR2 groups, where the haustoria did not successfully penetrate into host root and establish normal vascular continuity ( Figure 8B). In consideration of the fact that R proteins generally caused hypersensitive response (HR) characterized by programmed cell death (PCD), we asked whether they induce HR in transient expression assay. When these genes were overexpressed on Nicotiana benthamiana leaves via agro-infiltration, overexpression of CdNLR1 and CdNLR2 did not cause any visible HR on the infected area ( Figure 8C). Furthermore, DAB staining and measurement of cell death index also indicate that CdNLR1 and CdNLR2 expression did not induce HR in leaf tissue ( Figures 8D and  8E). Subcellular localization analysis by GFP-fused expression of CdNLR1 and CdNLR2 in an N. benthamiana transient expression system showed that they localized to the cell membrane ( Figure 8F). We speculated that CdNLR1 and CdNLR2 caused HR only in root but not in leaf tissue.
To further investigate the role of CdNLR1 and CdNLR2, we made transgenic N. benthamiana overexpressing these two genes underCaMV 35S promoter, namely CdNLR1-Ox and CdNLR2-Ox plants respectively. Both genomic DNA PCR with 35S promoter specific primers (upper lanes) and RT-PCR with CdNLR1 or CdNLR2 specific primers (middle lanes) confirmed the corresponding transgenic lines (Figures 9A and 9B). As O. cumana could parasitize on Solanaceae species including N. benthamiana, six-leaf stage transgenic plant seedlings were parasitized with O. cumana, and 50 days after parasitization O. cumana was removed from the soil for further analysis on parasite development. The number of emerging above-ground parasites on CdNLR1-Ox and CdNLR2-Ox plants decreased by 64% and 38%, and the underground parasites decreased by 69% and 44% respectively, compared with those on WT (wild type) plant ( Figures 9C-9F). We should also note that overexpression of CdNLR1 and CdNLR2 did not cause any abnormal phenotype, especially HR, in

DISCUSSION
Cross-species RNA movement was reported between host plant and the interacting organisms. Artificial RNA transfer from host to pathogen was intensely reported for many pathogenic organisms such as nematodes, insects and fungi in HIGS (host-induced gene silencing) experiment. 27 Recently, host mediated iScience Article RNAi toward pathogen gene was also verified in parasitic plants. 11,28,29 The above literature demonstrated the effective movement of ectopically expressed siRNA, a hallmark of RNAi, in host-pathogen interactions. Compared with this, natural RNA transfer event was rarely reported and the functional relevance of them were even less investigated. If any, such studies were mainly reported for small RNAs such as nat-siRNA (natural siRNA) and miRNA. 30 Accordingly, it was shown recently that sRNAs moved via vascular structures between host and pathogens including bacteria, fungi and parasitic plants. 12,13,31,32 Several pathogens deliver sRNAs into host cells to inhibit the immunity. For example, a set of 22 nt miRNAs could be delivered from the parasitic plant Cuscuta campestris to the hosts A. thaliana and N. benthamiana to suppress expression of target genes TIR1, AFB2, AFB3, BIK1 and SEOR1 to establish successful parasitism. 12 Furthermore, large-scale sRNA exchanges were also reported for Arabidopsis-Cuscuta parasite system. 33 Compared with sRNAs, it is difficult to investigate mRNA transfer due to the technical hurdles for distinction of long transcript sequence between source and destination organisms. 2,34 In this regard, parasitic plants offer a good option for the study of mobile mRNAs as they are evolutionarily distant from their host, and therefore sufficient sequence variation allow them to be distinguished between host and parasite. In early attempt, mRNA transfer was firstly confirmed for limited genes from host to parasite unidirection. 14 Largescale identification of mobile mRNAs was achieved via microarray analysis and sequencing technology for model hosts Arabidopsis, tomato and well-studied parasite dodder. 14,15 It is difficult to distinguish the mRNA between source and destination tissues, especially for non-model species which generally lack genomic data 35 (Westwood & Kim, 2017). However, now this issue is beginning to be overcome with the increases of sequencing depth and the improvement in bioinformatic approaches. Using multiple sequence alignment and comparative analysis, we identified around ten thousand mobile mRNAs between woody plant H. ammodendron and root-attached C. deserticola. In Arabidopsis-cuscuta connection, mobile mRNAs of 9,518 and 8,655 were originated from Arabidopsis and cuscuta respectively, accounting for 45% and 24% of the transcripts from source organism. 15 However, in the same report, only 347 (1.6%) and 288 (1.1%) mobile unigenes were identified in tomato-cuscuta parasitic system. 15 More recently, Liu et al. reported that 172 (0.06%) and 1,416 (0.28%) mobile mRNAs were exchanged between Arabidopsis and dodder, while mobile mRNAs were only 64 (0.12%) and 708 (0.19%) in soybean-dodder connection. 24 We report here that 2,370 (3.38%) and 7,496 (9.66%) unigenes moved between H. ammodendron and C. deserticola, respectively. Therefore, the scope and transfer specificity of RNA trafficking varied between different parasite systems. The present study represents a first report of mRNA transfer for root parasite and it differs from stem parasitic plant dodder. More confluent transfer of mRNA in CD/HA direction than the opposite direction might be contributable to species differences, however, it needs further investigation to be confirmed.
The most challenging issue for mobile RNA study is experimental verification of transfer event and functional characterization of mobile RNAs. Transcriptome sequencing is in depth and large scale, and confirmation of transfer mRNA could be facilitated by multiple sequence alignment as the host and parasite are distantly related to each other. 15,19 Here, we showed that sequence comparison with the orthologs from closely related species aided to identify the origin of mobile mRNA, whereas ruling out the sequence from distantly related species which represent the destination organism. Although the biological function of mobile mRNA is well studied for variety of genes within an organism, their implications in plant-parasite connections are poorly iScience Article understood. 2 This is mainly hampered by the lack of effective research tool for mobile mRNA, and is especially true for the system in which transgenic approach is lacking. To date, functional verification of mobile mRNAs was only achieved in model parasite system such as Arabidopsis-dodder and soybean-dodder in which transgenic methods are available, and for well-studied genes such as herbicide phosphinothricin acetyltransferase (PAT) 36 and florigen FT gene. 37 Here, we used sunflower-O. cumana system to facilitate both experimental validation and functional study of mobile mRNAs and confirmed their roles in parasite development. Overexpression of four C. deserticola genes CdAPM, CdATPA, CdPGI and CdRPL22 was shown to expedite haustorium attachment, indicating that these four genes possessed promotive roles in parasite attachment and haustorium development. Although less studied, the involvement of above genes for haustorium development was marginally reported. For example, several V-type ATP synthase/ATPase proteins were enriched in extrahaustorial space of powdery mildew-infected barley leaves to enhance disease susceptibility. FgGPI, a Glucose-6-phosphate isomerase gene which is involved in ATP biosynthesis and carbon source utilization, is required for hyphal growth and fungal development of Fusarium graminearum. 38 Taken together, our study demonstrates the effectiveness of sunflower-O. cumana parasitic system for the functional validation of C. deserticola genes and this strategy could also be useful for other root parasites.
It should be noted that C. deserticola is a specialist parasite that has only one reported host H. ammodendron, which is a sand-grown shrub of limited biotope in nature. 17 From an evolutionary aspect, excessive parasitism of C. deserticola would cause growth failure of host and renders both host and parasite to unfavorable growth situation. Therefore, it is reasonable to speculate that this narrow host-parasite relationship in natural ecosystem might evolve an equilibrium mechanism by which the parasite on a given host is restricted to a limited number for successful fulfillment of life cycle, instead of killing the host for their own death. Of interest, overexpression of CdNLR1 and CdNLR2, two C. deserticola NBS-LRR genes encoding putative R proteins, caused dramatic decrease in tubercle development. Therefore, transfer of PCDinducing mobile mRNAs for CdNLR2 and CdNLR3 might determine individual parasite destiny and prevents their over-parasitism on a host. Although the hypersensitive response was seldom reported for host-parasitic plant interactions compared with other host-pathogen system, it might also contribute to host immunity against parasitic plants. In stem parasite Cuscuta spp. and tomato association, a hypersensitive-like response (HLR) was reported to occur for incompatibility and to be mediated via plant hormone salicylic acid (SA

Limitations of the study
In this study, we identified around ten thousand mobile mRNAs between a woody host and a root parasite. However, it remains to be determined whether all these mobile mRNAs have functional implications in host-parasite interactions. Furthermore, lack of genetic manipulation for such non-model species hinders functional verification of the mobile mRNAs in original plant species; rather, we verified their roles in heterologous system. Therefore, the present study provides some evidence for cross-species mRNA transfer event whereas the individual mRNA role await further study. The functional correlation of transfer event and phenotypic outcome in recipient species will be unveiled in future study.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following: This study did not generate new unique reagents.

Data and code availability
The Illumina and PacBio raw sequence data for C. deserticola and H. ammodendron reported in this study have been deposited in the NCBI Sequence Read Archive (SRA) underBioproject PRJNA768605 (https:// www.ncbi.nlm.nih.gov/bioproject/PRJNA768605). Other existing, publicly available data used in this paper can be seen in method details. This paper does not report original code.
Data reported in this paper and any additional information required to reanalyze the data will be provided from the lead contact upon request.
The data that support the finding of this study are available in the supplementary material of this article.

Plant materials
Cistanche deserticola (Y. C. Ma) and Haloxylon ammodendron (C. A. Mey) were grown in a desert field at Dengkou, Bayinnor (40 45' N, 106 24' E) in China. The region has a typical continental arid climate with an average annual precipitation of 182 mm, about 70% of which occurs from July to September. Mean annual temperature of this region is 8.5 C, and mean temperature of the coldest (January) and the hottest (July) months are À15.7 C and 30 C, respectively. C. deserticola seed was evenly spread around and on top of the root of host H. ammodendron seedlings. At 80 days post inoculation, the parasite achieves the elongation stage and were used for transcriptome sequencing.
Sunflower cultivar (HN3638) were grown in growth chambers, and the Orobanche cumana (race G)-sunflower parasitic system was established in our lab. 25 Briefly, O. cumana race G seeds were collected in bulk from a highly infested field in Hohhot (China; latitude, 41.09289 N, longitude, 111.45785 E) during autumn 2018. The same sunflower cultivar (HN3638) that was heavily parasitized with O. cumana in the field was used for inoculation as it is highly susceptible to the O. cumana material we used here. Sunflower seed was sown in soil in pots and grown for 20 d to reach a height of $8 cm. The soil was carefully removed to expose host root, and surface-sterilized O. cumana seeds (100 mg per host) were evenly spread around and on top of the root tissue.

Bacterial growth
For the Agrobacterium-mediated transient expression, the cells of Agrobacterium tumefaciens strain GV3101 harboring the appropriate plasmids were cultured in LB broth containing appropriate antibiotics for vector selection, and the broth cultures were grown at 28 C and 200 rpm. See Method details for details.

Sample collection, library construction and transcriptome sequencing
The fresh succulent stems of C. deserticola (abbreviated as CD hereafter, and the same for other samples) 1 cm away from the haustorial junction was collected. At the same time, the roots of non-parasitized H. ammodendron (HA), the roots of CD-infested H. ammodendron (HC), and the haustorial interface (HI) were collected. Tissues were carefully harvested to avoid any chance of mixture and were rinsed for 10 min in 70% ethanol. Total RNA was extracted by Trizol reagent from individual plant tissue. RNA concentration was analyzed by NanoDrop2000 spectrophotometry and the RNA integrity value (RIN) was determined by Bioanalyzer 2100 system. Three biological replicates were performed for CD, HA and HC, and one replicate was used for HI. Ten mg of the total RNA, pooled from five individuals, was used for library construction, and RNA-seq was performed according to the methods described in our previous publications. 55

Transcriptome assembly and annotation
Raw reads in FASTQ input files were processed through in-house Perl script (https://www. perlscriptsjavascripts.com/perl/ba_faq/index.html). Adapter sequence, ploy-N reads and other low-quality reads were removed from raw data to obtain clean reads. Q20, GC content and sequence duplication iScience Article transfer unigenes were further separated with dual BLASTN (E value = 1e À10 ) against a collection of public available datasets from Orobanchaceae and Chenopodiaceae species. For Chenopodiaceae, we used the following datasets: H. ammodendron, next-generation transcriptome GSE63970 and GSE93684, 43,44 fulllength non-reduntant transcriptome in this manuscript; EST and mRNA sequences downloaded from NCBI database. For Orobanchaceae, we used the following databases: transcriptome of Triphysaria versicolor, Striga hermonthica and Phelipancheaegyptiaca, 57 EST and mRNA sequences downloaded from NCBI database. The putative CD/HC transfer mRNAs were compared with Orobanchaceae datasets and the low similarity sequences were removed. Similarly, the putative HC/CD transfer mRNAs were compared with Chenopodiaceae datasets, and the low similarity sequences were removed. At the same time, the highly conserved orthologs between HC and CD which also had high similarity with both Chenopodiaceae and Orobanchaceae were filtered out.

PCR confirmation of mobile mRNAs
Total RNA was isolated by Trizol reagent and genomic DNA contamination was removed by DNase treatment. RNA quality and concentration were detected by NanoDrop2000 spectrophotometer. GoScript reverse transcription system (Promega, A5000) was used for cDNA synthesis. Genomic DNA was isolated and RNA contamination was removed by RNase treatment. Both cDNA and genomic DNA were used as the templates for PCR amplification. The primers used were shown in Table S1. PCR reaction was performed, and the amplicons were separated on 1% agarose gel, with DNA marker run on the same gel for the estimation of band size. Representative amplicons from genomic DNA and cDNA were sequenced to verify the sequences.

pBI121 vector construction, transformation and parasite infestation
The ORF sequences for mobile mRNA genes were amplified from CD cDNA with the gene specific primers (Table S1). The amplicons were first inserted into cloning vector and were verified by sequencing. The sequences were then subcloned into pBI121 vector. GFP (Green fluorescent protein) sequence was amplified from pJL24 vector and C-terminally fused to these genes. Recombinant pBI121 vectors were transformed into competent cell of Agrobaterium tumefaciens strain GV3101 respectively, and cultured cells (OD 600 = 0.5) were suspended in the inoculation solution (10 mM MES, 10 mM MgCl 2 , 200 mM AS, 5% sucrose). Inoculation to sunflower seed was carried out according to the seed soak agroinfiltration method established in our lab. 25 Twenty daysold seedlings were parasitized with 100 mg of sterilized O. cumana seeds which were evenly spread around and on top of the root of sunflower seedlings. The seedlings were cultured at the growth chamber, and the parasites were taken from the soil for RT-qPCR assay, histological characterization and phenotypic observation. For subcellular localization, ORF sequence for full length CdNLR1 and CdNLR2 was amplified and was ligated into pEarleyGate103-SL vector. The recombinant vector was transformed to Agrobacterium tumefaciens GV3101, infiltrated into N. benthamiana leaves, and GFP signal (excitation at 488 nm) was visualized at 36-48 hpi. DAPI-stained nuclei were imaged using excitation at 405 nm and fluorescence detection at 425-493 nm.

QUANTIFICATION AND STATISTICAL ANALYSIS
All statistical analysis in this paper were analyzed by one way ANOVA between different groups using GraphPad Prism 8.